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Abstract 



A density-functional approach is used to calculate the inhomogeneous vortex 
density distribution in the flux liquid phase at the planar surface of a layered 
superconductor, where the external magnetic field is perpendicular to the su- 
perconducting layers and parallel to the surface. The interactions with image 
vortices are treated within a mean field approximation as a functional of the 
vortex density. Near the freezing transition strong vortex density fluctuations 
are found to persist far into the bulk liquid. We also calculate the height of 
the Bean-Livingston surface barrier. 
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Fluctuations of vortices in high-T c superconductors give rise to a number of novel effects 
IHIQ] the most prominent of which is probably the melting of the vortex lattice into a vortex 
liquid which was first proposed by Nelson || in 1988. While much work has been done to 
study the bulk properties of the vortex system, the influence of fluctuations on the behavior 
of flux vortices at the surface of a superconductor only rarely has been taken into account. 
Recently Brandt, Mints, and Snapiro 0] calculated the interaction of a single stack of ther- 
mally fluctuating pan-cake vortices with the planar surface of a layered superconductor. 
They found that fluctuations shift the short-range attraction between the vortices and their 
mirror-images to a long-range dipole-dipole interaction. In this Letter we consider a related 
problem, namely the interaction of a vortex-liquid with a planar surface. 

In order to describe the vortex system in the liquid regime we will apply the density- 
functional theory of classical liquids [El. The theory for the bulk vortex liquid has been 
worked out by Sengupta et al. |J who have calculated bulk correlation functions and treated 
the freezing of the vortex liquid into the solid phase. In the present work we will compute 
the inhomogeneous density distribution of the vortices in front of a planar surface using the 
Wertheim-Lovett-Mou-Buff ( WLMB) equation |7|] together with a mean-field treatment of 
the image-forces. [|J The same method has recently successfully been applied in a different 
context, namely the structure of liquid water at a metallic interface ||. From the equilibrium 
vortex-density profiles at the surface we calculate the local magnetic field which together 
with the image interactions determines the potential energy barrier JlO| for pan-cake vortices 
entering into the sample. From this the location of the irreversibility line (IL) in the B- 
T-diagram is estimated for the case of zero bulk pinning and vanishing geometrical surface 



barrier Wl 



Before discussing the surface problem we will shortly reconsider the density functional 
approach for the bulk vortex liquid following Sengupta et al. ||. This enables us to compute 
the bulk direct correlation functions we need as input into the surface calculations described 
below. In the limit of zero interlayer Josephson coupling and for H || c the vortices in an 
iV-layer superconductor can be viewed as a two dimensional classical system of N species 
of point particles interacting via the (asymptotically) logarithmic pair potential V n (r) of 
two vortices separated by n layers. In the limit iV — > oo the Fourier transform of this pair 
potential is given by |T2| 



e A 2 k\ + (Aid 2 ) sin 2 (M/2) 
V(k) = r L - 1 n n , (1) 



A 2 fc 2 +(4A 2 /d 2 )sin 2 (M/2) 

where k z (k±) is the component of the wave- vector k perpendicular (parallel) to the layers. 
The parameter eo = c?$q/47tA 2 is the characteristic energy scale, (3 is the inverse temper- 
ature, d <C A is the distance of the layers, and $ the flux quantum. The (2D-) pair 
potentials are repulsive for pancakes residing in the same layer and attractive but by the 
factor (d/2X) exp(—nd/X) weaker for pancakes in two separate layers of distance nd. How- 
ever, the attractive and repulsive logarithmic potentials compensate each other on distances 
larger than A leading to the usual short-range repulsion (~ K (r/X)) of Abrikosov vortices 
between two straight stacks of pancakes. 

From the pair potentials V n (r) one can derive the liquid-state pair correlation functions 
g n (r) within standard liquid theory. As in Ref. || we use the hypernetted chain (HNC) 
approximation. The basic equations are then the Ornstein-Zernicke equation 
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h n (r) = c n (r) + £?b W d 2 r'h m (r')c n „ m (\r - r'|) (2) 

and the HNC-closure relation 

g n (r) = exp [-(3V n {r) + h n (r) - c n (r)] (3) 

for the total correlation function h n (r) = g n (r) — 1 and the direct correlation function 
c n (r), where ^ is the bulk density of the vortex liquid. For the numerical (iterative) 
solution of these equations it is appropriate to reformulate Eq. (^) in k-space and to separate 
the direct correlation function c n (r) into a long-range and a short-range part. (Note that 
c n ij) — > —f3V n (r) as r — > oo.) We were able to obtain self consistent solutions of the 
complete set of Eqs. (0) and ([3]) for a periodic system with 512 layers using a Fast Hankel 
transform method [[R|. 

The phase boundary between the vortex liquid and the vortex solid was obtained in a 
different way than in Ref. || namely by examining the thermodynamical stability of the 
system, i.e. we looked at the spinodal curve rather than the coexistence line at the first 
order phase transition. The system becomes unstable if the second functional derivative of 
the grand canonical potential w.r.t. the density, A4(r, z, r', z') = 5 2 Q/5g(r, z)5g(r f , z'), is no 
longer positiv definite. In the case considered here Ai(r, z, r', z') is diagonal in k-space 



because of the translational invariance of the vortex-system in all three space directions, 
-M(k, k') = (1 — c(k z , k±))S(\c — k'). As expected, near the spinodal curve the system 
becomes unstable w.r.t. to fluctuations of the density ~ e* k±r with |kj_| « ^ng 1 ^ 2 (k z = 0), 
which indicates that the formation of a vortex lattice becomes favorable. The freezing line 
(spinodal curve) in the (£>,T)-diagram is then obtained by extrapolating 1 — c(k z ,k±) to 
zero. For high magnetic fields B the transition temperature approaches the 2D-melting 
temperature T 2 d which is independent of the vortex density . 

We now turn to the calculation of the local vortex density near the surface of the super- 
conductor. We start with the Wertheim-Lovett-Mou-Buff (WLMB) equation J7| which is an 
exact relation following from classical density-functional theory. In a more general context 
this system of integro-differential equations relates the particle densities g a of species a to 
an external potential V a : 

Vi log g a {l) = -/3ViK(l) + J2 I d2V 2 Qf){2)c a p{l, 2, [g]) (4) 

a 

Here, c Q( g(l,2, [g]) is the inhomogeneous direct correlation function which is a functional of 
the density itself. In the following we will again work in the HNC approximation and replace 
c aj g(l, 2, [g]) by its bulk version. We introduce a cartesean coordinate system where the (y,z)- 
plane is the surface of the superconductor which fills the half-space at x > 0. Because of the 
translational invariance in ^-direction we are left with one equation for the vortex density 
g(x) in each layer which is a function of x only. Now, the bulk direct correlation function 
c n (r) is written as a sum of a long-range and a short-range part, c n (r) = c£ R (r) — (3V n {r), 
we carry out the integration parallel to the surface, and define 

oo 

^R,\\(„\ _ f 2rdr SR , 
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With En V n (r) = (e /2ir)K (r/\) the WLMB-equation §) then becomes 



e -\x-x'\/\_ 



\ogg(x) = pAg\e e- x / x + PF\x,[g]) + / dx'^c 5ii 'll(|x - x'\) - ^f- f dx'-£ 

J ax 2 J ax 

o o 

(5) 



Eq. (^j) is the basis for our numerical calculation. The first term on the r.h.s. of Eq. (^) 
gives the force produced by an exponentially decaying shielding current at the surface of the 
superconductor which arises from the homogeneous solution of the London equation. The 
second term is the image force which is the sum of the self-image force F SI and the force 
poi p roc l ucec i by the images of all other pancakes, [g]) = F SI (x) + F OI (x, [g]). The 

third term in Eq. (||) gives the mean force caused by the liquid correlations and the last 
term represents the electromagnetic vortex-vortex interaction. 

Eq. (||) has to be solved for the density g(x) at x > with g(x) J — ? gs while the 
parameter Ag is determined by the condition that the potential of mean force, Vpmf(x) = 
— ksT ln(g(x) / gs) is zero at the surface (f?(0) = qb)- The local magnetic field B(x) then 
reads 

B(x) = $ (A£ - g B )e- x/x + g(x')e-\ x '-^ x dx' } (6) 

o 

the applied field is B a = B(0), and the bulk magnetic field is given by B = B(oo) = gs&o- 
In the following we will derive an expression for the image forces F 1 within a mean field 
approximation. The self image force acting on a pancake vortex at position r = (x, 0) is 
given by 

F SI (x) = -^ (x>0, (7) 

where the cut-off £ is related to the coherence length. The ^-component of the force exerted 
on the same vortex by the image of a second vortex located at r' in layer n is 



V 

/n(r,r') = 



n (\/ r i2 + 4xa/) (x + x') 



r? 9 + Axx' 



with r\2 = |r — r'|. Thus, the average total force produced by the non-self images can be 
written 

F°\x, [£]) = £/ / n (r, v') e {v')h n {x'\v)dx' (8) 

where h n (r'\r) is the total correlation function for the inhomogeneous liquid distribution at 
the surface. Here we assume that forces parallel to the surface cancel each other on the 
average. If the total correlation function h n (r'\r) is again approximated by its bulk version 
h n (ri2) the non-self image force F 01 can finally be expressed as a linear functional of the 
vortex density g(x), 
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F OI (x, [q]) = J Q(x')J(x,x')dx'. (9) 



The kernel 3{x, x') is given by 



_ , 2r 12 h n (r 12 )V; Ur{ 2 + Axx') 
J{x,x') = {x + x')J2 / V ' dr 12 . (10) 

" > 'r 2 2 — (x — x') 2 Jr 2 2 + 4xx' 



\x-x'\ 

It is important to notice that the force of the non-self images cancels the long-range self- 
image force for sufficient large x. Indeed, writing V n (x) ~ g n m(x/£) for x > A with 'charges' 
go = —1, Qn^o ~ (d/2\)e~ nd / x (Y^nQn = 0), setting ^(a;) = and using the 'screening 
condition' for the total correlation function h n (r), q + J2nl d 2 rg B q n h n {r) = 0, one can 
perform the integrals (|9|) and ( |10"D analytically and gets F 01 — — ► eo/47rx = —F SI . 

In what follows we present results of the numerical solution of Eq. ([5]). We introduce 
the dimensionless parameters / = \ 2 Qb, T = /3eo, 5 = d/X, and n = A/£. All calculations 
have been performed for k = 100 and 5 = 0.01 while / and T were varied. About 10 4 
iterations of Eq. (^) on a grid of 2048 sites were needed to achieve convergence. However, 
when approaching the freezing transition of the vortex liquid the rate of convergence becomes 
slower and the numerical procedure finally diverges. The kernel J{x, x') for the image-forces 
(Eq. fllCf )) can be calculated once and then stored in a lookup-table. 

Fig. 1 shows the so-obtained vortex density profiles for / = 0.75 and parameters V = 100, 
250, and 360. Freezing occurs at approximately T = 380. It is seen that near the freezing 
transition strong density fluctuations persist far into the bulk liquid, so that the vortex-liquid 
is still solid-like at the surface. These fluctuations are washed out at higher temperatures 
(lower values of T). The corresponding magnetic field near the surface is shown in the left- 
hand inset of Fig. 1. Directly at the surface (invisible in Fig. 1) the density rapidly drops 
from to almost zero, which reflects the existence of a Bean-Livingston surface barrier 
ITJ . This surface barrier is shown in the right-hand inset of Fig. 1 where the potential of 



mean force Vj>mf is plotted for several values of the parameter / = 0.6,0.75,1.0,2.0 and 
T = 250. The barrier height U turns out to be nearly independent of the temperature but 
decreases with increasing magnetic field (parameter /). This is shown in the inset of Fig. 2. 

When the applied field is increased above its equilibrium value vortices will start to enter 
into the sample, where the flux-entry rate 71 is determined by the Boltzmann factor e~^ u . 
The irreversibility line (IL) separates the irreversible and the reversible part of the phase 
diagram on experimental time scales and is thus given by TZ = const. We can therefore 
estimate (up to a constant factor) the location of the surface-barrier induced IL in the 
.B-T-diagram by the condition U(B) = aksT, where a is a constant. 

In order to compare with experimental results we will in the following use typical pa- 
rameters for Bi 2 Sr 2 CaCu208 (BSCCO). The London penetration depth is taken to have a 
temperature dependence A(T) = A(0)(1 - (T/T c ) 4 ) 1/2 with A(0) = 1400 A and T c = 93 
K, k = 100, and d — 15 A. Fig. 2 (broken line) shows the melting line (spinodal curve) 
of the flux-line lattice in BSCCO as obtained by the method described above. For T > 50 
K it agrees well with the experimental data of Majer et al. [|rT| (o) . Deviations at low 
temperatures are possibly an effect of the Josephson interaction which has not been taken 
into account here. Fig. 2 (solid line) also shows the estimated location of the IL. A constant 
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factor a = 42 was chosen such that the curve lies near the experimentally obtained IL (+) 
from Ref. [If for a sample of triangular cross-section, where the geometrical surface barrier 



is absent. Since bulk pinning plays no role in these experiments either, as argued in Ref. ITT 



the IL shown in Fig. 2 originates from a Bean-Livingston surface barrier alone. Though our 
simple calculation may not be applicable to this experiment in detail because of the more 
complicated sample geometry and possible effects of the surface roughness, the qualitative 
features of the IL found here are the same as those found in the experiment. In both cases 
the IL has a negative slope dB/dT and is much steeper than the melting line crossing the 
latter one at some point. 

In conclusion we have calculated the density distribution of pan-cake vortices in the 
flux liquid phase at the surface of a layered superconductor using density functional theory 
and a mean field treatment of the image-interactions. Near the melting line strong density 
fluctuations are found to reach far into the bulk liquid. This could possibly enhance bulk- 
pinning near the surface. It is shown that a Bean-Livingston barrier exists from which we 
estimate the location of the irreversibility line (IL) in the S-T-plane which is qualitatively 
consistent with experimental findings for BSCCO crystals. 

One of us (A.K.) thanks the Deutsche Forschungsgemeinschaft for financial support. 
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FIGURES 



FIG. 1. Density profile g(x) for constant bulk density (/ = 0.75) and different temperatures 
(r = 100, 250, and 360). The right-hand inset shows the potential of mean force (PMF) at constant 
temperature (r = 250) for different bulk densities (/ = 0.6,0.75,1.0,2.0, from the upper to the 
lower curve). The left-hand inset gives the internal magnetic field B{x) at the surface for the same 
set of parameters. 

FIG. 2. The melting line of the vortex lattice in Eh^S^CaC^Os: Experimental data from 
Ref. [11] (o) and this work (dashed line). Irreversibility line (IL): Experimental data for a sample 
with triangular cross-section from Ref. [11] (+) and as estimated in this work (solid line, see text). 
The inset shows the barrier height U/eq as function of /. 
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